2.2 PPAC 信号处理-II¶
1. Tracking¶
利用多个PPAC的x，y信息进行traking的原理如图所示,图中红点代表探测器测量的位置。由于PPAC的探测效率，对于每个入射粒子只有部分PPAC的x或y方向可给出位置信息。假设对于某一个入射粒子，有两个以上(含两个)探测器给出x位置信息 (如图:1Ax,1Bx,3x) 时，入射粒子的x-z径迹可由上述x-z点进行线性拟合得到。y-z径迹的处理方法与x-z径迹的处理方法一致。利用上述线性方程，可内插或外推某一z位置的x和y信息。图中空心黑点代表由线性方程计算得到的位置。
本实验中采用的重建策略为：¶
1) 若F8PPAC3中有有效的位置而且F8PPAC1和F8PPAC2中至少有一层PPAC有有效的位置。
2) 若F8PPAC3没有有效的位置，则
- a) 这个事件若被 F8PPAC1和F8PPAC2中的3层或3层以上的PPAC探测。
- b) 若一个事件只有2层PPAC探测到，要求这两层PPAC分别来自两套PPAC1和PPAC2。
束流径迹-利用所有PPAC的信息¶
2. PPAC 的位置分辨率¶
准直束刻度方法¶
通过准直孔将粒子准直，入射到探测器上。此时测得的位置分辨率$\sigma$与探测器本征位置分辨率$\sigma_0$的关系为 $$ \sigma^2=\sigma^2_0+\sigma^2_{geo} $$ ，其中$\sigma_{geo}$为准直孔引入的分辨率（图a）。
上述方法一般用放射源进行，而PPAC探测器的位置分辨与入射粒子的种类以及能量相关，因此由放射源得到的分辨率不能直接用在束流条件下。
- 束流条件下，如果在探测器前放置准直孔，则由于束流在准直孔上产生散射(图b)，导致$\sigma_{geo}$，不仅与孔的几何条件有关，也与束流的条件相关。
多探测器径迹测量方法(束流)：¶
- 选择被全部PPAC(5层)测到的事件进行径迹重建，得到每个探测器的残差$\Delta x^i=x^i-x^i_{track}$的分布，其中$x^i$和$x^i_{track}$分别为测量和拟合得到的位置。径迹外推到靶所在的位置，得到靶点坐标$tx_{track}$。
- 每个探测器的残差分布的宽度$\sigma^i_{\Delta x}$是所有探测器本征分辨率的函数，即$\sigma^i_{\Delta x}=\sigma^i_{\Delta x}(\sigma^0_0,\sigma^1_0...,\sigma^4_0)$。
- 探测器的本征分辨率可通过Monte Carlo模拟得到。具体做法如下：对于一组给定分辨率$(\sigma^0_0,\sigma^1_0...,\sigma^4_0)$，每个探测器的位置$x'^{i}$ 由高斯抽样 $Gauss(x^i_{track},\sigma^i_0)$给出。拟合径迹得到每个探测器的拟合位置$x'^i_{track}$,以及残差$\Delta x'^i=x'^i-x'^i_{track}$分布，径迹外推到靶点的坐标为$t_{x'}$。调节不同的分辨率,使得模拟和实验的残差分布符合。这样就能间接算出每一层PPAC的位置分辨率。此时靶点位置残差$\Delta t_{x'}=t_{x'}-t_{x'track}$分布的宽度即为束流在靶上的位置分辨率。
3. PPAC 探测效率¶
假设有$N$个粒子穿过探测器灵敏面积，探测器实际探测到的数目为$N_{det}$, 则探测效率$\epsilon$为 $$ \epsilon = \frac{N_{det}}{N} $$
在本实验设置下，可以有很多方法求探测器的效率：
例1.ppac的阳极信号的计数作为分母，阴极计数作为分子（图a）。
例2.选择任意两组(i,j)探测器，拟合得到径迹, 记录径迹穿过探测器k灵敏面积的所有事件数目$N_{ij}$。记录在上述条件下探测器k具有正确位置信息的事件数目$N_{kij}$（图b）。 $$ \epsilon = \frac{N_{kij}}{N_{ij}} $$
用上述做法可分别计算探测器的x，y，x-y的效率
示例代码¶
下面示例代码假设同时有1A，2A，3探测器给出x，y的位置信息。
Branch¶
$PPACF8[i][j]$ -Calibrated Data(位置已刻度好)
- x,y的信号有效时，同时要求anode信号有效。
- 不正确的位置信息用-999表示。
| Branch | PPAC | code |
|---|---|---|
| PPACF8[0][0] | PPAC 1 Layer A X (mm) | xx[0] |
| PPACF8[0][1] | PPAC 1 Layer A Y (mm) | yy[0] |
| PPACF8[0][2] | PPAC 1 Layer A Z for X-plan (mm) | xz[0] |
| PPACF8[0][3] | PPAC 1 Layer A Z for Y-plan (mm) | yz[0] |
| PPACF8[0][4] | PPAC 1 Layer A Anode time (ns) | |
| PPACF8[1][0] | PPAC 1 Layer B X (mm) | |
| PPACF8[1][1] | PPAC 1 Layer B Y (mm) | |
| PPACF8[1][2] | PPAC 1 Layer B Z for X-plan (mm) | |
| PPACF8[1][3] | PPAC 1 Layer B Z for Y-plan (mm) | |
| PPACF8[1][4] | PPAC 1 Layer B Anode time (ns) | |
| PPACF8[2][0-4] | PPAC 2 Layer A * | xx[1],yy[1],xz[1],yz[1] |
| PPACF8[3][0-4] | PPAC 2 Layer B * | xx2b,yy2b,xz2b,yz2b |
| PPACF8[4][0-4] | PPAC 3 * | xx[2],yy[2],xz[2],yz[2] |
利用makeclass生成 tracking.h, tracking.C¶
root -l f8ppac001.root
>tree->MakeClass("tracking");
>.q
编辑 tracking.h和tracking.C 保存
使用方法：
root -l
> .L tracking.C
> tracking t
>t.Loop()
>.q
root -l tracking.root //分析
tracking.h¶
在traking.h 代码中增加了用户变量，成员函数SetBranch()，TrackInit() 以及SetTrace()定义
#ifndef tracking_h
#define tracking_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
// Header file for the classes stored in the TTree if any.
class tracking {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
Int_t fCurrent; //!current Tree number in a TChain
// Fixed size dimensions of array or collections stored in the TTree if any.
// Declaration of leaf types
Float_t PPACF8[5][5];
Float_t F8PPACRawData[5][5];
Int_t beamTrig;
Int_t must2Trig;
Float_t targetX,targetY;
//by user
Double_t xx[3],xz[3],yy[3],yz[3];//1A,2A,3
Double_t xx2b[2],yy2b[2],xz2b,yz2b;//2B x,y, 0-measured, 1- fitted.
Double_t anode2b;
Double_t dx[3],dy[3];//residual
Double_t tx,ty;//target position
Double_t c2nx,c2ny;//chi2/ndf for xfit,yfit
// List of branches
TBranch *b_PPACF8; //!
TBranch *b_F8PPACRawData; //!
TBranch *b_beamTrig; //!
TBranch *b_must2Trig; //!
TBranch *b_targetX; //!
TBranch *b_targetY; //!
tracking(TTree *tree=0);
virtual ~tracking();
virtual Int_t Cut(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry);
virtual Long64_t LoadTree(Long64_t entry);
virtual void Init(TTree *tree);
virtual void Loop();
virtual void SetBranch(TTree *tree);//by user
virtual void TrackInit();//by user
virtual void SetTrace(TH2D *h,Double_t k,Double_t b,Int_t min,Int_t max);
virtual Bool_t Notify();
virtual void Show(Long64_t entry = -1);
};
#endif
#ifdef tracking_cxx
tracking::tracking(TTree *tree) : fChain(0)
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
if (tree == 0) {
TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("f8ppac001.root");
if (!f || !f->IsOpen()) {
f = new TFile("f8ppac001.root");
}
f->GetObject("tree",tree);
}
Init(tree);
}
tracking::~tracking()
{
if (!fChain) return;
delete fChain->GetCurrentFile();
}
Int_t tracking::GetEntry(Long64_t entry)
{
// Read contents of entry.
if (!fChain) return 0;
return fChain->GetEntry(entry);
}
Long64_t tracking::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
if (!fChain) return -5;
Long64_t centry = fChain->LoadTree(entry);
if (centry < 0) return centry;
if (fChain->GetTreeNumber() != fCurrent) {
fCurrent = fChain->GetTreeNumber();
Notify();
}
return centry;
}
void tracking::Init(TTree *tree)
{
// The Init() function is called when the selector needs to initialize
// a new tree or chain. Typically here the branch addresses and branch
// pointers of the tree will be set.
// It is normally not necessary to make changes to the generated
// code, but the routine can be extended by the user if needed.
// Init() will be called many times when running on PROOF
// (once per file to be processed).
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fCurrent = -1;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("PPACF8", PPACF8, &b_PPACF8);
fChain->SetBranchAddress("F8PPACRawData", F8PPACRawData, &b_F8PPACRawData);
fChain->SetBranchAddress("beamTrig", &beamTrig, &b_beamTrig);
fChain->SetBranchAddress("must2Trig", &must2Trig, &b_must2Trig);
fChain->SetBranchAddress("targetX",&targetX,&b_targetX);
fChain->SetBranchAddress("targetY",&targetY,&b_targetY);
Notify();
}
Bool_t tracking::Notify()
{
// The Notify() function is called when a new file is opened. This
// can be either for a new TTree in a TChain or when when a new TTree
// is started when using PROOF. It is normally not necessary to make changes
// to the generated code, but the routine can be extended by the
// user if needed. The return value is currently not used.
return kTRUE;
}
void tracking::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
if (!fChain) return;
fChain->Show(entry);
}
Int_t tracking::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns 1 if entry is accepted.
// returns -1 otherwise.
return 1;
}
#endif // #ifdef tracking_cxx
tracking.C¶
在tracking的基础上进行了修改
#define tracking_cxx
#include "tracking.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TFitResult.h>
void tracking::SetBranch(TTree *tree)
{
//measured pos
tree->Branch("xx",&xx,"xx[3]/D");//1A,2A,3
tree->Branch("xz",&xz,"xz[3]/D");
tree->Branch("yy",&yy,"yy[3]/D");
tree->Branch("yz",&yz,"yz[3]/D");
//difference between measured and calculated -for pos resolution.
tree->Branch("dx",&dx,"dx[3]/D");
tree->Branch("dy",&dy,"dy[3]/D");
//2B x,y,anode
tree->Branch("xx2b",&xx2b,"xx2b[2]/D");
tree->Branch("yy2b",&yy2b,"yy2b[2]/D");
tree->Branch("anode2b",&anode2b,"anode2b/D");
//target x-y
tree->Branch("tx",&tx,"tx/D");
tree->Branch("ty",&ty,"ty/D");
//ch2/ndf for linear fitting.
tree->Branch("c2nx",&c2nx,"c2nx/D");
tree->Branch("c2ny",&c2ny,"c2ny/D");
tree->Branch("beamTrig",&beamTrig,"beamTrig/I");
tree->Branch("must2Trig",&must2Trig,"must2Trig/I");
tree->Branch("targetX",&targetX,"targetX");
tree->Branch("targetY",&targetY,"targetY");
}
void tracking::TrackInit()
{
tx=-999;
ty=-999;
//1A
xx[0]=PPACF8[0][0];
yy[0]=PPACF8[0][1];
xz[0]=PPACF8[0][2];
yz[0]=PPACF8[0][3];
//2A
xx[1]=PPACF8[2][0];
yy[1]=PPACF8[2][1];
xz[1]=PPACF8[2][2];
yz[1]=PPACF8[2][3];
//3
xx[2]=PPACF8[4][0];
yy[2]=PPACF8[4][1];
xz[2]=PPACF8[4][2];
yz[2]=PPACF8[4][3];
//2B
xx2b[0]=PPACF8[3][0];
yy2b[0]=PPACF8[3][1];
xz2b=PPACF8[3][2];
yz2b=PPACF8[3][3];
anode2b=PPACF8[3][4];
xx2b[1]=-1000;
yy2b[1]=-1000;
}
void tracking::SetTrace(TH2D *h,Double_t k,Double_t b,Int_t min,Int_t max){
if(h==0) return;
if(min>=max) return;
for(int i=min;i<max;i++){
h->Fill(i,(Int_t)(i*k+b));
}
}
void tracking::Loop()
{
TH2D *htf8xz=new TH2D("htf8xz","xz trace by ppac",2200,-2000,200,300,-150,150);
TH2D* htf8yz=new TH2D("htf8yz","yz trace by ppac",2200,-2000,200,300,-150,150);
TFile *opf=new TFile("tracking.root","recreate");
TTree *tree=new TTree("tree","ppac traking");
SetBranch(tree);
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
TrackInit();
bool b1a=abs(xx[0])<150 && abs(yy[0])<150;
bool b2a=abs(xx[1])<150 && abs(yy[1])<150;
bool b3=abs(xx[2])<150 && abs(yy[2])<150;
if(!b1a || !b2a || !b3) continue;
//fit x-z trajectory
TFitResultPtr r;
TGraph *grx=new TGraph(3,xz,xx);
TF1 *fx=new TF1("fx","pol1",-2000,0);
r=grx->Fit(fx,"SQ");
xx2b[1]=fx->Eval(xz2b);
tx=fx->Eval(0);
SetTrace(htf8xz,fx->GetParameter(1),fx->GetParameter(0),-1800,0);
for(int i=0;i<3;i++) dx[i]=xx[i]-fx->Eval(xz[i]);
c2nx=r->Chi2()/r->Ndf();
delete grx;
delete fx;
//fit y-z trajectory
TGraph *gry=new TGraph(3,yz,yy);
TF1 *fy=new TF1("fy","pol1",-2000,0);
r=gry->Fit(fy,"SQ");
yy2b[1]=fy->Eval(yz2b);
ty=fy->Eval(0);
SetTrace(htf8yz,fy->GetParameter(1),fy->GetParameter(0),-1800,0);
for(int i=0;i<3;i++) dy[i]=yy[i]-fy->Eval(yz[i]);
c2ny=r->Chi2()/r->Ndf(); //对于任何程序的自动拟合，原则上都要输出拟合误差进行评估
delete gry;
delete fy;
tree->Fill();
if(jentry%10000==0) cout<<"processing "<<jentry<<endl;
}
htf8xz->Write();
htf8yz->Write();
tree->Write();
opf->Close();
}
In [1]:
%jsroot on
TFile *ipf=new TFile("tracking.root");
TTree *tree=(TTree*) ipf->Get("tree");
TCanvas *c1=new TCanvas("c1","c1");
束流径迹¶
In [2]:
TH2D *hxz=(TH2D*) ipf->Get("htf8xz");
hxz->Draw("colz");
c1->Draw();
In [3]:
TH2D *hyz=(TH2D*) ipf->Get("htf8yz");
hyz->Draw("colz");
c1->Draw();
束流在靶上投影¶
- 假设靶与束流线垂直。
In [4]:
tree->Draw("ty:tx>>htx(120,-60,60,120,-60,60)","must2Trig","colz");
c1->Draw();
In [5]:
tree->Draw("tx:ty>>htx(120,-60,60,120,-60,60)","beamTrig","colz");
c1->Draw();
$\chi^2/Ndf : tx$
残差，chi2/NDF¶
残差分布并不是单一的gaus分布，有一个窄的gauss重叠在宽的gauss之上。 其可能的原因为：
- Nucl.Instr.Meth.A 593,376 (2008)
- The position dependent resolution is responsible for the multiple-Gaussian shape of the residual distributions presented in the preceding section. Consequently,the detector performance is expressed in terms of an optimal resolution corresponding to the narrow Gaussian residual width and an effective resolution associated with the area-weighted average of the two Gaussian widths.
In [6]:
TF1 *total2=new TF1("total2","gaus(0)+gaus(3)");
In [7]:
TF1 *g1 = new TF1("g1","gaus");
TF1 *g2 = new TF1("g2","gaus");
TF1 *total=new TF1("total","gaus(0)+gaus(3)");
TH1F *hdx;
Double_t sigma;
In [8]:
tree->Draw("dx[0]>>hdx(200,-5,5)");
hdx=(TH1F*)gROOT->FindObject("hdx");
hdx->Fit("g1","","",-0.5,0.5);
sigma=g1->GetParameter(2);
total->SetParameter(2,sigma);
total->SetParameter(5,5*sigma);//初始化，估计半宽为2*sigma；
hdx->Fit("total");
gPad->SetLogy();
c1->Draw();//residual
In [9]:
gPad->SetLogy(0);
tree->Draw("dx[0]:dx[1]");
c1->Draw();
In [10]:
tree->Draw("c2ny:dy[0]>>hh(40,-10,10,200,0,1000)","","colz");
c1->Draw();//从chi2/ndf图上可看出，部分事件的径迹拟合误差很大，这部分要在后续数据处理中去掉。
In [11]:
tree->Draw("c2ny>>hh(200,0,1000)","","");
gPad->SetLogy();
c1->Draw();//从chi2/ndf图上可看出，部分事件的径迹拟合误差很大，这部分要在后续数据处理中去掉。
In [12]:
gPad->SetLogy(0);
tree->Draw("tx:ty>>(120,-60,60)","c2nx<10 && c2ny<10 && beamTrig ","colz");
c1->Draw();//
In [13]:
tree->Draw("tx:ty>>(120,-60,60,120,-60,60)","(c2nx>20 || c2ny>20) && beamTrig ","colz");
c1->Draw();//
PPAC2B x,y,x-y的探测效率¶
In [14]:
Long64_t N_track;
TCut c2btrack="abs(xx2b[1])<100 && abs(yy2b[1])<100";//径迹穿过探测器的灵敏面积
Long64_t Na_det,Nxa_det, Nya_det,Nxya_det;
TCut c2banode="anode2b>-1000";//阳极有信号
TCut c2bxa="abs(xx2b[0])<100";//x面有正确信号，即x面和阳极同时有信号
TCut c2bya="abs(yy2b[0])<100";//y面有正确信号，即y面和阳极同时有信号
1. 计算穿过探测器灵敏面积的粒子数目¶
In [15]:
tree->Draw("yy2b[1]:xx2b[1]>>(200,-100,100,200,-100,100)",c2btrack,"colz");//fitted
N_track=tree->GetEntries(c2btrack);//得到给定条件下的计数。
cout<<N_track<<endl;
c1->Draw();
2. 计算在上述条件下，探测器的有效信号数目¶
In [16]:
tree->Draw("yy2b[0]:xx2b[0]>>(200,-100,100,200,-100,100)",c2bxa&&c2bya&&c2btrack,"colz");
Na_det=tree->GetEntries(c2banode && c2btrack);
Nxa_det=tree->GetEntries(c2bxa && c2btrack);
Nya_det=tree->GetEntries(c2bya && c2btrack);
Nxya2_det=tree->GetEntries(c2bxa && c2bya && c2btrack);
cout<<Na_det<<" "<<Nxa_det<<" "<<Nya_det<<" "<<Nxya2_det<<endl;
c1->Draw();
3. 计算anode,x, y, x-y 的探测效率¶
In [17]:
Double_t effa,effxa,effya,effxya2,effxya;
In [18]:
effa=Double_t(Na_det)/N_track;
effxa=Double_t(Nxa_det)/N_track;//ex*ea
effya=Double_t(Nya_det)/N_track;//ey*ea
effxya2=Double_t(Nxya2_det)/N_track;//ex*ey*ea*ea
effxya=Double_t(Nxya2_det)/N_track /effa;
In [19]:
TString eff;
eff.Form("PPAC2B:\n eff_a=%2.f%%,\n eff_xa=%.2f%%, \n eff_ya=%.2f%%",effa*100,effxa*100,effxa*effya*100);
cout<<eff.Data()<<endl;
eff.Form("\n eff_xya2=%.2f%% \n eff_xya=%.2f%%",effxya2*100,effxya*100);
cout<<eff.Data()<<endl;
In [20]:
!jupyter nbconvert 2.2_PPAC_tracking.ipynb --to html
In [ ]: